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Abstract 



We study the response of a large array of coupled nonlinear oscillators to parametric excita- 
tion, motivated by the growing interest in the nonlinear dynamics of microelectromechan- 
ical and nanoelectromechanical systems (MEMS and NEMS). Using a multiscale analysis, 
we derive an amplitude equation that captures the slow dynamics of the coupled oscil- 
lators just above the onset of parametric oscillations. The amplitude equation that we 
derive here from first principles contains uncommon nonlinear gradient terms which yield 
a unique wave-number dependent bifurcation similar in character to the behavior known 
to exist in fluids undergoing the Faraday wave instability. We suggest a number of exper- 
iments with nanomechanical or micromechanical resonators to test the predictions of our 
theory, in particular the strong hysteretic dependence on the drive amplitude. 
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Introduction 



In the last decade we have witnessed exciting technological advances in the fabrication 
and control of microelectromechanical and nanoelectromechanical systems (MEMS and 
NEMS). Such systems are being developed for a host of nanotechnological applications, as 
well as for basic research in the mesoscopic physics of phonons, and the general study of 
the behavior of mechanical degrees of freedom at the interface between the quantum and 
the classical worlds [1-3]. Surprisingly, NEMS have also opened up a new experimental 
window into the study of the nonlinear dynamics of discrete systems with many degrees 
of freedom. A combination of three properties of NEMS resonators has led to this unique 
experimental opportunity. First and most important is the experimental observation that 
micro- and nanomechanical resonators tend to behave nonlinearly at very modest ampli- 
tudes. This nonlinear behavior has not only been observed experimentally [4-11], but has 
already been exploited to achieve mechanical signal amplification and mechanical noise 
squeezing [12,13] in single resonators. Second is the fact that at their dimensions, the nor- 
mal frequencies of nanomechanical resonators are extremely high — recently exceeding the 
IGHz mark [14] — facilitating the design of ultra-fast mechanical devices, and making the 
waiting times for unwanted transients bearable on experimental time scales. Third is the 
technological ability to fabricate large arrays of MEMS and NEMS resonators whose col- 
lective response might be useful for signal enhancement and noise reduction, as well as for 
sophisticated mechanical signal processing applications. Such arrays have already exhibited 
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interesting nonlinear dynamics ranging from the formation of extended patterns [15] — as 
one commonly observes in analogous continuous systems such as Faraday waves — to that 
of intrinsically localized modes [16]. Thus, nanomechanical resonator arrays are perfect for 
testing dynamical theories of discrete nonlinear systems with many degrees of freedom. At 
the same time, the theoretical understanding of such systems may prove useful for future 
nanotechnological applications. 

This work is motivated by a recent experiment of Buks and Roukes [15], who succeeded 
in fabricating, exciting, and measuring the response to parametric excitation of an array of 
67 micromechanical resonating gold beams. Lifshitz and Cross [17] described the response 
of the beams with a set of coupled nonlinear differential equations ()2.4jl . They used secular 
perturbation theory to convert these equations into a set of coupled nonlinear algebraic 
equations for the normal mode amplitudes of the system, enabling them to obtain exact 
results for small arrays but only a qualitative understanding of the dynamics of large arrays. 
In order to obtain analytical results for large arrays we study here the same system of 
equations, approaching it from the continuous limit of infinitely-many degrees of freedom. 
Our central result is a scaled amplitude equation ()5.15|] . governed by a single control 
parameter, that captures the slow dynamics of the coupled oscillators just above the onset 
of parametric oscillations. This amplitude equation includes uncommon nonlinear gradient 
terms and exhibits a unique wave-number dependent bifurcation that can be traced back 
to the parameters of the equations of motion of the system ()2.4|) . We confirm this behavior 
numerically and make suggestions for testing it experimentally. 

This thesis is organized as follows. Chapters 1-3 serve as a brief background, summa- 
rizing the results of Buks and Roukes (chapter 1) and Lifshitz and Cross (chapters 2,3) 
that are essential for the understanding of the work presented in the following chapters. 
In chapter 4 we present a detailed description of our derivation of the amplitude equations 
describing the response of large arrays. A reduction to a single amplitude equation just 
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above the onset of oscillations is preformed in chapter 5. Single mode solutions of the 
amplitude equation are discussed in chapter 6, and experimental schemes for obtaining 
such solutions and their unique properties are suggested and demonstrated numerically in 
chapter 7. 
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Chapter 1 



Experimental Motivation 



This work is motivated by a recent experiment of Buks and Roukes [15, henceforth BR] who 
constructed an array of micromechanical resonating beams forming a diffraction grating. 
Fig- Hm shows a micrograph of 67 gold beams that BR fabricated on top of a silicon nitride 
membrane. Each beam was 270 x 1 x 0.25 fim in size, and the distance between two 
neighboring beams was 4/xm. The membrane was then removed, leaving only the ends of 
the beams connected to the silicon surface. Electrostatic forces between the beams were 
applied by connecting them alternately to two base electrodes. Such forces induce tunable 
coupling between the beams, and thus a collective spectrum of vibrational modes could 
be obtained. By introducing an ac component to the electrostatic forces, these modes can 
be parametrically excited. Optical diffraction was used to study the response of the array 
as a function of the driving frequency and the dc component of the electrostatic forces, 
Vdc- BR excited the array near its second instability tongue, i.e. the driving frequency lied 
within the band of normal frequencies (see chapter Ej). 

Before applying the electrostatic forces, the characteristics of each beam separately were 
measured. The averaged fundamental frequency was 179.3 kHz, with a standard deviation 
of 0.53 kHz, and the quality factors Q ranged from 2, 000 to 10, 000. No correlations 
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Figure 1.1: A side view micrograph of the array constructed by Buks and Roukes. The 
purple background is the silicon surface, from which a rectangular piece was etched away 
(black rectangular) , leaving the fabricated gold beams connected to the silicon surface from 
both edges. A potential difference V is applied between the two base electrodes. 

were found between the location of the beams in the array and their specific properties. 
Therefore in the following model (chapter [SJ the beams are regarded as identical. Once 
the electrostatic forces were applied, the quality factor values severely decreased as Vdc 
was increased. This suggests that the dissipation of the system is mainly due to induced 
currents between the beams. 

Fig. 11.21 shows the measured collective response of the array as a function of the dc 
component of the electrostatic forces and the driving frequency. For each value of Vdc the 
amplitude of the ac component was set to 50 mV and the frequency was gradually increased. 
The band of frequencies the array responded to became wider as Vdc was increased. This 
can be understood already at the linear analysis level, it is a direct consequence of the 
formation of the band of collective modes. The larger the dc component is, the wider 
is the band. The upper bound of the band however strongly depends on Vdc, while the 
linear analysis predicts it to be simply the fundamental frequency of an individual beam 
Ub/2TT ^ 179.3 kHz. Moreover, for certain values of the dc component the array responded 
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Figure 1.2: A color map showing relative intensities of the diffracted light from the array, 
as a function of the voltage Vdc and the driving frequency / taken from [15]. The dashed 
white line shows a fit to the measured lower bound frequency obtained from a simple 
linear theory. For Vdc ~ 14V^ the array responded at frequencies beyond the expected 
upper bound frequency ~ 179.3kHz. 



to frequencies beyond this upper bound frequency. It is also evident that as the driving 
frequency was swept up, the response showed a small number of wide peaks, much wider 
than the 67 resonance peaks that a linear theory predicts. These features were qualitatively 
explained by Lifshitz and Cross [17] by taking into account the nonlinearities of the array. 
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Chapter 2 



Model and Equations of Motion 

Lifshitz and Cross [17, henceforth LC] showed that the nontrivial response of the array is a 
direct consequence of the nonlinear character of the beams. They derived a set of equations 
of motion, introducing only the essential terms for capturing the nonlinear features obtained 
by BR. We shall briefly discuss some of their main guidelines for deriving the equations of 
motion. 

The normal frequencies of the individual beams were found to be well separated. Thus 
for moderate driving each beam is strictly vibrating in its fundamental frequency ui,. As- 
suming the beams are restricted to in-plane vibrations only, each beam can be described 
by a single degree of freedom Un, its displacement from equilibrium. 

For simplicity, the attractive forces due to the applied electrostatic potential between 
the beams are approximated by nearest neighbor interactions 

F^?Lic = -lrnu;lA^[l + H cos{2u;pt)]{un+i - 2u„ + (2.1) 

and A^H are the dc and ac components of the electrostatic forces respectively, where 
the prefactor mLofj2 is inserted for convenience. The parametric excitation introduced into 
the system by the ac component, is an instability of the system that occurs for driving 
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frequencies 2uj.p close to one of the special values luj^jn [18], where ujq is one of the normal 
frequencies of the array and n is an integer that labels the so-called instability tongues of 
the system. In the BR experiment, the system was excited in its second instability tongue, 
i.e. n = 1 ,ijjp ^ ojq/2. Since the response at the second tongue was found to be quite 
similar to the response at the first tongue, the calculations preformed by LC were done for 
the first instability tongue. In the current work we proceed with concentrating on the first 
tongue only, and thus we always take Up itself to lie within the normal frequency band. 

The restoring forces of the beams are due to their elasticity. Measurements done on 
individual beams indicate that the linear restoring force is supplemented with a cubic 
term of the displacement acting to stiffen the beam. Neglecting higher order nonlinear 
corrections the elastic force can therefore be written as 

f'^Lc = -rn^l («n + ««n) , (2.2) 

with a > 0. Motivated by the observations of BR regarding the quality factors (see 
chapter HJ LC assumed that the dissipation of the system is mainly due to the electrostatic 
interaction, which induces currents through the beams. Therefore the dissipation should 
depend on the difference variable m„ — m^-i, describing the relative displacements of a 
pair of neighboring beams. LC showed that a nonlinear dissipation term must also be 
introduced in order to obtain bounded response for driving frequency sweeps. Thus the 
dissipation terms were taken to be ^ 

Fdils = ^rnuJbT{un+i - 2m„ + tt„„i) 

+ ^mUJbarj[{Un+l - U„)^(m„+i - tt„) - {Un - Un-lf{Un - (2.3) 

^The effect of introducing gradient-dependent dissipation terms instead of local dissipation terms is to 
renormalize the bare dissipation coefHcients F and rj to wavenumber-dependent coefficients of the form 
4rsin^ {q/2) and 4r/sin^ (9/2) respectively, as will become apparent below. 
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The dissipation of the system is assumed to be weak, which makes it possible to excite 
the beams with relatively small driving amplitudes. In such case the response of the 
beams is moderate, justifying the description of the system with nonlinearities up to cubic 
terms only. The weak dissipation can be parameterized by introducing a small expansion 
parameter e <^ 1, physically defined by the linear dissipation coefficient T = ej, with 7 of 
order one. The driving amplitude is then expressed by A^H — eh, with h of order one. 
The weakly nonlinear regime is studied by expanding the displacements Un in powers of 
e. Taking the leading term to be of the order e^/^ ensures that all the corrections, to a 
simple set of equations describing coupled harmonic oscillators, enter the equations at 
the same order of e^/^. 

Introducing the scaled variables t t/ujh and Un — > Un/ y/a, LC eventually obtained 
the following set of dimensionless equations 

Un + Un + ul+^ [A^ + eh COs{20Jpt)] {Un+l - 2Un + Un-l) 

- ^T] [{Un+l - UnYiUn+l " O " (^n " Un-lf {Un " Un-l)] = 0. (2.4) 

The boundary conditions are set according to the experiment of BR, who had two additional 
fixed beams at both ends of the array, thus uo — un+i — 0. 
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Chapter 3 



Normal Modes and Linear Stability 



Analysis 



We first expand the response of the beams into N standing waves, 

N 

Un = y~] ^mjt) sin(g^n), with Qm = j;^—^, m = l...N. (3.1) 

m=l 

Substituting ()3.H) into the equations of motion ()2.4|) yields N nonhnear ODEs for the 
amplitudes of the waves, 



where N.L. stands for nonlinear terms that couple the N linear Mathieu [19] equations, 
and LJm is given by the dispersion relation 



+ + e2 sin^(g„/2) 7$^ - h cos(2cjpt)<l>„ + N.L. = 



(3.2) 




1 -2A2 sm\qj2). 



(3.3) 
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One can easily verify that the zero-displacement state, m„ = for all n, is always a solution 
of the equations of motion (j2.4jl . though it is not always a stable one. To study the 
transition to oscillating solutions, we first linearize the equations of motion about the zero- 
displacement state, and introduce a new time scale T = et upon which the growth of small 
perturbations will occur. T and t are regarded as independent variables, following the 
multiple scales procedure [18], as we expand the response in orders of e, 

$^(t,T) = ei/2$^„(t,T) + e3/2$„^(t,T) + 0(e^/2), (3.4) 

and express the time derivatives in terms of both t and T, 

dt^dt + edr and dl dl + 2ed^T + ^(e^). (3.5) 

Substituting (|3.5p into the linear part of Eq. ()3.2|) and collecting terms of same order e we 
obtain for the e^^^ order terms 

{dl + ujl)^^, = Q, (3.6a) 

and for the e^/^ order terms 

{dl + ujI) = [-2d^T - 27 sin^ (9-/2) dt + h sin^ (g„/2) (e^^-.* + ^-^2^.t^-^ _ 

(3.6b) 

It is convenient to write the solution of Eq. ()3.6a|) in the following form 

$,^0 = A{T)e"'-' + C.C. , (3.7) 
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where A{T) is a complex amplitude which can slowly change over long time scales. Sub- 
stituting into Eq. ()3.6hjl we obtain 

{dl + ul) = - (220;^^ + 2iiJml sin^ {qj2) A^^ e^"-* (3.8) 
+ h sin^ A^e*(2^*'+^™)* + h sin^ (g„/2) A*^e'^^'^^-'^^'>' + c.c. 

The terms on the right hand side of Eq. (j3.8j) proportional to e*'^™* (called secular terms) 
excite the equation for $^1 in its resonant frequency. In order to prevent the perturbative 
term from diverging, we must account for a solvability condition, demanding that the 
sum of all the secular terms vanishes^. For large frequency detunings Up — Um ^ e, only the 
first two terms of the left hand side of Eq. ()3.8p must vanish and the solvability condition 
calls 

dA 

^ = -7sin2(g„/2)A„, (3.9) 

yielding a slow exponential decay of Am and thus no excitation of the m*'* mode. On the 
other hand for frequency detunings of order e the fourth term of the right hand side of 
Eq. ()3.8|) is also secular. We introduce a detuning parameter cUp = Um + |e^m, and the 
solvability condition yields 

- (2tiUm^ + 2tuJm7 sm^{qm/2)Am] + h sm''{qj2)Al/^"^^ = 0. (3.10) 



dT 

Trying a solution of the form Am{T) = ame°"'"^e*^^, with a™ G K yields 

-2iUmCrmam + UJm^mam " 2iuJm1 SlYl^ {qm/'2)am + h SVD? {q^ / 2) a*^ = 

{um^rnf + (2u;„7 sin^ (g^/2) + 2a;^(T„)^ = h"^ sin^(g„/2) 



cxm = -7sin^(gm/2) ± \nhsiYi\qm/2) /2umf - (f^™/2)2. (3.11) 



^The constraint on the secular terms can be obtained in a wider context of linear differential operator 
theorems, as explained briefly in section ITU . 
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If the linear growth rate am > the solution will grow, thus the critical driving amplitude 
hm{uJp) is the one for which the growth rate am vanishes, explicitly, 

hmiuJp) = 2Uml\ 1 + .^7 = "^^mlj 1 + (^^^^^^r^) ■ (3-12) 

y \2-fsm\qm/2)J y \e-f sm\qm/2) J 

hm{ujp) is the driving amplitude required to excite the m*^ normal mode, from zero dis- 
placement into parametric oscillations with frequency Up. Note that it has the familiar 
form of a first instability tongue [19], modified by the dispersive correction sin^(gm/2). 

3.1 The Case of Distinct Normal Frequencies - A Single 
Mode Response 

If the spacing between the normal frequencies 6uJm = ~ ^m-i is much greater than e 
one should expect the system to respond like a single degree of freedom. Only the mode 
for which Up is close to its normal frequency is excited. Fig. 13.11 (A) shows five hm{ojp) 
curves in the (/i, Up) plane which mark the onset of standing waves with a wave number 
qm for an array of 5 resonating beams. The curves do not overlap, indicating that in this 
array for moderate drives only single modes are excited. 

In order to calculate the steady state response, a single mode solution should be sub- 
stituted into the equations of motion (I2.4jl 

Un = ^m{t) sm{qmn). (3.13) 

By taking into account the nonlinear terms neglected up to this point, the saturation of 
Am can be calculated. LC performed such a calculation and obtained the steady state of 
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the single mode solution 



Un = e^^'^2\am\cos{ujpt — ip)sm{q.mn), with 

(3.14a) 

1 /g \2 2 

fl^ = TTT I -|«mP - t^mf^m I + {2Urnl + QuJ^Tj ii\V? {q„J 2)\arn\'^) ,and 



1 /2sin^(gm/2)7 + 6sin'^(gm/2)?7|amp 
V9 = -arctan 9, |2_,, o 



(3.14b) 



(3.14c) 



In Fig. 13.11 (B) and (C) the response of the excited m*'^ mode la^p is plotted as a function 
of the amplitude h for two values of the frequency detuning Qm- Solid curves indicate 
stable solutions and dashed curves unstable solutions. The nature of the response changes 
significantly as the detuning is set above a critical value of flc = 16 sin^(gm/2)u;m7?7/3. 
For Qm < the amplitude of the oscillating solution grows continuously for h above 
threshold h > hm{ojp) and is stable. This is known as a supercritical Hopf bifurcation [20]. 
On the other hand, a subcritical Hopf bifurcation is obtained for Vt^ > ^c- An unstable 
solution grows below threshold, until it reaches the so-called saddle-node point, where 
the curve of I function of h bends around (a turning point), and the solution 

becomes stable and an increasing function of h. The stability of the steady state solution 
is determined by linearizing the solutions (j3.14jl about a small perturbation with the same 
spatial dependence sin{qmn) (see section lOl bellow) . Whether these perturbations decay 
or grow determines the stability of the solution. The stability can alternatively be deduced 
from a corollary of the factorization theorem [21], which states that the stability of the 
solutions switches either at turning points or at points at which two solutions intersect. In 
both approaches only instabilities towards the growth of perturbations with the same wave 
number Qm are examined, thus such a stability analysis is valid only for cases of sufficiently 
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separated normal frequencies. 

3.2 The Case Overlapping Stability Curves 

If is sufficiently large such that 5uJm — ^^jj ~ ~ several modes can be excited 
simultaneously for moderate driving amplitudes h ~ 0(1). This situation is illustrated for 
= 100 oscillators in Fig. 13.21 which shows in (A) a set of overlapping instability tongues 
hmiyjp) for 64 < m < 81 plotted as a function of t^p, and in (B) the values of hmiyjp) 
as a function of mode number m for a particular value of Up = cjys (dashed line in (A)), 
outlining the neutral stability curve below which the zero-displacement state is stable. 

In the limit A^ — > cxo of very large arrays, the frequency spectrum becomes essentially 
continuous and so does the neutral stability curve of Fig. 13.21 (B). The first mode qc to 
emerge when increasing the driving amplitude from zero will then be the mode which 
minimizes the continuous hmiyjp) curve, occurring at a critical driving amplitude he- For 
a moderate-size system, with a discrete frequency spectrum, the first excited mode will 
be the one whose critical driving amplitude hm{ujp) is closest to h^. By further increasing 
the driving amplitude, the zero-displacement state becomes unstable to a band of wave 
numbers bounded by the neutral stability curve, as indicated by the horizontal line in 
Fig. 13.21 (B). Once the excited modes start saturating into standing waves they interact 
with each other through the nonlinear terms, potentially yielding complicated dynamical 
behavior, as observed by LC [17] in the exact solutions obtained using secular perturbation 
theory. The exact solutions however could only be obtained for small arrays, enabling only 
a qualitative understanding of the dynamics of large arrays. 
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Figure 3.1: (A) The five hmiojp) curves of an array of 5 beams in the {h,ujp) plane are 
shown. Inside the region bounded by the rn}^ curve the zero displacement state is unstable 
towards the growth of the rn}^ mode. The dashed vertical lines mark two excitations of the 
array with frequency detunings of = 0.1 (blue line) and = 2 (red line). The critical 
detuning is f2c = 0.196. (B) For < f2c a supercritical bifurcation is obtained , (C) while 
1^4 > VLc yields a subcritical bifurcation. The parameters used are A = 0.4, r/ = 0.1 and 
e7 = 0.001. 
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Figure 3.2: (A) A set of overlapping instability tongues hmij^p) for an array of 100 oscil- 
lators as function of ti^p. Modes 64 < m < 81 are plotted, the dashed vertical line indicates 
the value of Up used in B. (B) Neutral stability curve as a function of mode number m, 
plotted for ujp = ^^73. Marks the minimal driving amplitude hmi^^p) required to excite 
the m'^^ mode from the zero-displacement state. The horizontal line indicates the band of 
unstable modes for a particular value of h = 2. The parameters used are A = 0.5 and 
e7 = 0.01. 
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Chapter 4 

The Response of Large Arrays - The 
Amphtude Equation Approach 

In order to study arrays with a large number of beams, we adopt an approach usually 
applied to continuous spatially-extended systems. In such systems the emergence of modes, 
accessible from the unstable band, appear as modulations of the basic pattern determined 
by the initial unstable mode qc- These modulations develop slowly over long spatial scales 
and thus can be described by an envelope function of the basic pattern whose dynamics 
obeys an appropriate set of amplitude equations [22]. 

4.1 Derivation of The Amplitude Equations 

The general form of the amplitude equations can be deduced from symmetry consider- 
ations, but here we wish to derive the equations explicitly and obtain exact expressions 
for all the coefficients, that can subsequently be tested quantitatively both numerically 
and experimentally. We introduce a continuous displacement field u{x,t), keeping in mind 
that only for integral values n of the spatial coordinate does it actually correspond to the 
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displacements u{x = n,t) = Unit) of the discrete set of oscillators in the array. We also 
introduce slow spatial and temporal scales, X = ex and T = et, upon which the dynamics 
of the envelope function occurs, and expand the displacement field in terms of e, 

u = e^/2M(o)(x, t, X, T) + €''^u^^\x, t, X, T) + 0{e'/^). (4.1) 

In order to substitute the expansion ()4.ip into the equations of motion ()2.4p we must 
express the temporal and spatial operators in Eq. p.4|l in terms of the new variables T 
and X. The time derivative substitution is given by Eq. (i;i5|) . The displacement of the 
nil beam is given by 

nn±i{t) = e'/'u^'\x±l,t,X,T) + Oie'/'). (4.2) 

We insert Eqs. ()3.5|] . (14.111 and ()4.2|1 into ()2.4p . and collect terms of the same order of e. 
The e^/^ terms yield 

£«(°) = 0, £ = dl + ^A^ (e^- + e"^- - 2) + 1, (4.3) 

where e^"^ is an operator that performs a spatial shift by 1, thus the full spatial operator 
in Eq. (j4.3j) is simply the discrete Laplacian. Eq. (I4.3jl is solved by setting 

u^°\x,t,X,T) = {A+{X,T)e-''^^'' + A*_{X,T)e''^^^) e''^^' + c.c. (4.4) 

The response to the lowest order of e is therefore expressed in terms of two counter- 
propagating waves with modulated complex amplitudes and A_. This is a typical 
ansatz for parametrically excited systems, though the choice of the scaling of X and T 
with e is not. The asterisk and c.c. stand for the complex conjugate and Qp is determined 
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by the driving frequency Up and the dispersion relation ()3.3p 



'1 — uj"^ 

gp = 2arcsiny (4.5) 

In order to obtain the e^/^ order contribution to the different terms in the equations of 
motion, we express the 0(e^/^) correction of Eq. ()4.2|1 in terms of u^^^ and A± using the 
substitution 

A^{X + e, T) ^ A±(X, T) + e^^^^l^. (4.6) 
Using Eqs. ()4.2|) . ()4.4|) and ()4.6p these contributions are 



(8 A dA* \ 

Un+l + - 2m„ : 



U 



(1) 



(x + 1, t) + u^^\x - 1, t) - 2M(^)(a;, t) 



'(9/4 BA* 
aX aX 



u 



(1) 



(a; + 1, t) + M^^^a; - 1, i) - 2uW(x, t) 



9 A* \ 

- 2i sin(g„) — ie"*^''^' - ^e'^^"" e*"*'* + c.c. ; (4.7b) 

\ oX oX J 

e(^n+i + Un-i - 2u„) : 

- Aiup sin2(gp/2) (A+e'^^"^ + A*_e*''''^) e*""* + c.c. ; (4.7c) 



ecos(2u;pt)(M„+i + Un~i - 2un) ■■ 
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2sin2(gp/2) (A+e^^«''^ + A^e*^^^) e'^'^^^'^^^* 

2 sm2(gp/2) (^;e'«''^ + A_e-'«''^) e^^^^^"^^)* + c.c. ; (4.7d) 



+ O (e^^'"^*, e^^'^^^) + C.C. ; (4.7e) 

-— ((A+e-^''^'^(e^^*^' - 1) + Ale'«''^(e^^«^ - l))e^'^''* + c.c.)^ 
dt 

X (A+e-'«''^(e^**^ - 1) + A:^e^*^^(e^^*^ - 1)) e^"^ + 0(e^^'^''*, e^^«^^) + c.c. 
= [ \e-''i^ - l\^{\A+\^ + \A_\^) + yl+^_e-'29p^(e^*'''' - 1)^ 

+ ^;^^e^2«^^(e^*«^ - 1)2 ] (A+e-^*^^(e^*«^ - 1) + Ale^«^^(e^*«^ - 1)) e^^^ 

+ 0(e^3'"^*,e^3''''^) + C.C. ; (4.7f) 

and 

- Unfiiln+l - iln) + (^n-l " Mn)^(M„-l - : 

iup\e'''^^ - l|2(e^«'' + e-''^^ - 2){\A+\^ + 2|>l_|2)>l+e-^*^^e^'^''* 
+ icUp|e-^«^ - l|2(e*«^ + e-''^^ - 2)(2|^+|2 + \A_\^)A*_e''i^''e'^^^ 
^ 0(^(^^vt^ eisqpx^^ + c.c. 

= -tl6u;pSm\qp/2) [ {\A+\^ + 2|^_n^+e-''^-"e*"''*+ (4.7g) 
{2\A+\^ + |A_|2)Ale*«''^e^'"^*] + 0(e^=^'"^*, e^=^«^^) + c.c. , 

where (^(^e*^'^''*, e'^*^) are terms proportional to e'^'^''* or e*^*^ which do not enter the 
dynamics at the lowest order of the e expansion. The e^/^ order terms of the equations of 
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motion ()2.4|1 therefore yield 



^U^'^ = f 

+ [h {Ale''^"'' + A.e"*''^^) - i2up-f (A+e"*'''^ + A*_e''"''')] sm\qp/2)e''^^' 

- (3 + iSriiUp sm\qp/2)) {\A+\^ + 2\A^\^) A+e-'^^^ e''^^' 

- {2> + iSr]UpSin\qp/2)) (2 1 A+l^ + 2) A^e^^^^e*"^* + 0(6^=^^^*, e^^'?''^) + c.c. 

(4.8) 

Let us denote by Vq the zero eigenvectors of the self adjoint operator £, namely, e~*(''p^'~'^p*) 
and e^iip^+'^p^) and their complex conjugates. If m*^^^ is a solution of Eq. ()4.8p . then / must 
be orthogonal to the zero eigenvectors because 

(/, Vo) = V,) = {u^^\ ZV,) = 0, (4.9) 

where (/,(?) oc / f{x,t)g*{x,t)dxdt is the inner product of the L2{C) space. The Fredholm 
alternative theorem [23] states that Eq. (j4.9ll is also a sufficient condition for the existence 
of a solution u^^^ for Eq. ()4.8p . Thus the solvability condition of Eq. ()4.8|) requires that the 
coefficients of the zero eigenvector terms at its right hand side must vanish identically. We 
therefore obtain two coupled amplitude equations, 

^ + "^9^ = -l^^^\qp/m+ - sin2(g,/2)A_ (4.10a) 

- (^4r/sin^(g,/2)- z^^ (|A+p + 2|A_p) 
dA_ dA^ h 

'W~'''~dX ^ -7sin'(W2)A_ +z— sin2(gp/2)A+ (4.10b) 

- (^r^s\n\qp/2)+i^^ (2|A+p + l^.^) 
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where Vg = ^ = 2^ is the usual group velocity. Eqs. M.lOjl are two coupled complex 

Ginzburg-Landau equations (CGLE) [24], similar to the amplitude equations previously 
derived for describing Faraday waves excitations [25,26]. 



4.2 Linear Stability Analysis of The Zero-Displacement 
State 

We now wish to reexamine the stability of the zero-displacement state, using the amplitude 
equations (|4.1()ll we derived. For small perturbations of the zero displacement state 



A± = a±(T)e' 



-ikX 



\a±\ <^ 1, 



(4.11) 



only the linear terms of Eqs. (I4.1()jl are considered. Thus we can write the amplitude 
equations in a matrix form 



d_ 



( \ 













ivgk --fsin^{qp/2) -i^sm'^{qp/2) 
\ ^2^sin2(gp/2) -Wgk--fsm^{qp/2) ^ 



( \ 

« 1 



(4.12) 



describing the dynamics at the onset of modes with a wave number q = qp + ek. The 
solutions of Eq. (I4.12|) can be expressed as linear combinations of two exponents in time, 
Be"+'^ and De"-^ , with growth rates a± determined by the eigenvalues of the matrix in 
the right hand side of the equation. These are given by 



cr±{h) = -7sin^(gp/2) ± J {hsin\qp/2) /2upY - {vgkf. 



(4.13) 



23 



The critical amplitude hk{ujp) required to excite oscillations with a wave number qp + ek is 
obtained by setting the larger eigenvector cr+(/ifc) = 0, thus 



The minimum of hk{ujp), he = '^^ojp, is obtained for A; = 0, i.e. the amplitude equations 
approach yields an initial instability with a wave number qc = qp. Comparing Eq. ()4.14|) 
with the neutral stability curve Eq. (I3.12|l obtained by a direct linear stability analysis of 
the equations of motion ()2.4|1 reveals the scope of accuracy of the scalings we use. 

For <^ e, the spectrum of the normal modes can be regarded as being essentially 
continuous, by replacing and uJm by continuous variables q and uj related by a continuous 
dispersion relation ()3.3p . In this limit, with Up taken from the normal frequency band, q^ 
is determined by the minimum of the continuous neutral stability curve Eq. ()3.12|1 , 



Even though the beams oscillate at frequency cjp, the initial unstable wave number qc is 
only approximately qp. We focus on small detunings Vt = 2{ujp — uj)/e <^ 1 and expand 




(4.14) 




(4.15) 



Eq. sum in n, 





(4.16) 



1^^ 2ujp^ 2\2jsm\q/2) 



( en 1 f n 




+ o{n\en^). 
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Identifying the frequency detuning with 



= ujp-uj = \vg\{q - qp) = e\vg\k, 



(4.17) 



and comparing Eq. ()4.16|1 with the expansion of Eq. ()4.14p in small k reveals that in the 
amplitude equation approach the wave number detuning ek corresponds to q — qp up to 
order corrections. This suggests that the scaling we chose for the spatial coordinate 
X = ex is valid for the emergence of waves with g — ^ e^. 

For h = hkiujp) the two eigenvalues are cr+ = and cr_ = — 27sin^gp/2, with the 
corresponding eigenvectors being 



I 1 ^ 



and F_ 



I 1 ^ 



where e 



7 sin^(gp/2) — ivgk 



\/^^sm\qp/2) + {vgky' 

(4.18) 

The amplitudes B and D can be expressed in terms of a± by multiplying Eq. ()4.12|) by 
the matrix 

/ \ ' / 

111 i I e-'"^ -1 

-e'^ 1 



J 



2 s\.n{ip) 

from the left. We then obtain that B and D are given by 



V 



(4.19) 




(4.20) 



and obey 



d_ 
df 



-27sm=fe/2) 



(4.21) 



as expected. At the onset of oscillations with a wave number qp + ek the linear combination 
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B becomes unstable while D decays exponentially with a growth rate (7_ = —2^sm'^{qp/2), 
suggesting that the dynamics close to the onset of oscillations can be well described using 
a single complex amplitude. This is what we do in the next section. 
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Chapter 5 

Reduction to a Single Amplitude 
Equation 

Motivated by the linear analysis of Eqs. ()4.1()jl preformed in chapter 14 .21 we express the 
amplitudes A± in terms of B and D, 





\ 


_/ 1 


1 \ 




\ 


U- 


/ 






u 


/ 



Just above threshold h — he he the D amplitude decays exponentially in time. After 
transients that last over periods of order T^^ ~ e^^, the response can therefore be expressed 
by the B amplitude only 





\ 


_/ 1 


\ 


U- 


/ 




I 
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B{X,T) = \B{X,T)\e'^B{x,T) ^presents an envelope of a single standing wave as apparent 
from substituting ()5.2|] back into (j4.4|) 



u 



(o)(a;,t,X,T) 



{B{X, T)e-''^^'' + B{X, T)*e^(^''^-^)) e''^"* + c.c. 



2\B{X,T) \ (cos {upt - qpX + ips) + cos {upt + qpX - ips - i')) 



4:\B{X, T)| cos(c<;pt - ij/2) cos{qpX - ^b{X, T) - 



(5.3) 



^/2 is the temporal phase relative to the drive. 

This reduction of the description of the dynamics to a single amplitude B is similar to 
the procedure introduced by Riecke [27] for describing the onset of Faraday waves. It is 
typical of parametric driving, which excites a single standing wave at threshold. 

5.1 Scaling of A± Just Above Threshold 

We define a reduced driving amplitude g by letting 



In order to obtain an equation, describing the relevant slow dynamics of the new amplitude 
-B, we need to select the proper scaling of the original amplitudes A±, as well as their spatial 
and temporal variables, with respect to the new small parameter 5. 

If the coefficient of nonlinear dissipation i] is of order of one, it is apparent from the 
original amplitude equations M.lfljl that the cubic terms saturate the growth of the am- 
plitudes A±. However, it is physically realistic to assume that rj is small. Therefore a 
quintic term must enter in order to saturate the growth of the amplitudes A±. This can 
be achieved by defining the small parameter 6 with respect to the coefficient of nonlinear 



{h - he) /he = gS, 



5< 1. 



(5.4) 
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dissipation, letting 



V 



51/2 



r/o ~ 0{1] 



(5.5) 



and taking the amplitudes to be of order S^^'^. 

It is apparent from Eq. ^4. I'M for the linear growth rates that with a drive amplitude 
that scales like 6 the growth rate of the amplitude B scales like 6 as well. The bandwidth 
of unstable wave numbers scales as 5^/^, as obtained by Eq. ()4.14|] . The new temporal and 
spatial scales are therefore defined hy t = 6T and ^ = respectively, and we finally 

make the ansatz that 



(a: 




(A 








[' J 



^«;W(X,T,e,r) ^ 



v 



V 



«(X,T,e,' 



+ 55/4 



J 



(5.6) 



The phase e''^ defined by ()4.18|) to the lowest order of 5 is i, because the band of unstable 
wave numbers k scales like 5^/^. 



5.2 Derivation of The B Amplitude Equation 

The amplitude equation for B{^,t) is derived by once again using the multiple scales 
method. We substitute Eqs. ()5.4p . ()5.5|) and the ansatz (|5.6|) into the original amplitude 
equations ()4.in|) . and collect terms of the same order of S. The 5^^^ order terms satisfy 
Eqs. (I4.1()|) trivially, due to the fact that the 5^^'^ order terms in the ansatz ()5.6|] were 
chosen according to the linear stability analysis of Eqs. ()4.inp performed in section IT!2l 
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In order to obtain the higher order terms in S we perform the following calculations 



\A+\^A. 

\A^\^A 

dA+ _ 
'df~ ~ 
dA_ 



dT 

dA+ 
~dX 

dA^ 



+ = 6'/'\B\'B + 6'/\B'w^'^* + 2\B\V) ■ 

+ = 6'/'\B\'B + 6'/\iB%^'^* - + \B\'w^'^) ■ 

_ = 6^/'\B\'tB + 55/4(_52^,(i)* + 2|i?|\;«) ; 

.3/4^^ e5/4^^ , .5/4^ . 

^ dT ^ dr' 

3/4 9^ x3/4^ , .5/4^^ , .5/4^^ . 

dx ^ dx ^ 9e ' 



(5.7) 



Collecting all the 5^^^ order terms of Eqs. ()4.inp yields 



— z 



(5.8) 



where D is a linear operator given by 



D 



dr + Vgdx + 7 sin2(gp/2) 
-Z7 sin^(gp/2) 



z7sin^(gp/2) ^ 



dr - Vgdx + 7 sin^(gp/2) 



(5.9) 



There is no solvability condition since the right hand side of Eq. (I5.8jl is automatically 
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orthogonal to the zero eigenvalue of D, ^ ^ ^ . The solution of Eq. ()5.8p is given by 



w 



dB . 9 



-Vn— + i—\B\'B 



27sin^(gp/2) V 2cUp 



f 1 ^ 



(5.10) 



We plug Eq. ()5.1flp back into Eqs. ()4.inp using Eqs. ()5.7|) . collect all the terms with ^^Z'* 
and obtain 
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[Ot + Vgdx + 7sin2(gp/2)] w^^^ + i^i siTi^{qp/2)v^^^ 



dB 



^ +7sm2(gp/2)^S - 12r)osm\qp/2)\B\''B 



+i A ( B'w^^y + 2\Bfw'^'^ + 2iB^v^'^* - 2z|5|2^;W + 2\B\^w^'^) 
2ujp \ / 



dB 



d 



Or ^sm\qj2)di \ 2 4^;^ 



-Vg dB 



+ i—\B\'B] + 'ysm\qp/2)gB - 127josm\qp/2)\B\^B 



+i 



27 sin {qp/2)uj 



p I 



B' 



-Vg dB* 



. 9 



4|5p 



-Vn dB 



+ i- 



9 



2 (9{ Au. 



\B\^B 



+2iB'^ -i 



dB ^ 



.Vg dB* 
'2~dC 

vl 



+ 



AuJr. 



\B?B* 



.VadB 9 



d^B 



dr ^ 2-fsm^{qp/2) d^^ ULUp-fsm^{qp/2) 



-12rjosm\qp/2)\B\^B + 



2i\B\'' i^^ + —\B\^B 



9i;„ 



+ 7sin'(?p/2)^S 



u;p7sin {qp/2) 



-3 + 



Vn 



ujp-fS\v?{qp/2,) 



3 3\ „2^S 



dB 



4 ' 2; 9e a;27sin2(g,/2) V 8 2 



27 27 27 ^27^^ 
4 4 



2>Vn 



2ujp-fs\v?{qp/2) 



m 



,dB 



+ B' 



.dB* 



81 



9e ; 8a;2^sm2(gp/2) 



(5.11a) 
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and 



dB dv^^^ 

-i— + v„^T^ + sm'^{qp/2)giB ~ 12r]osm^{qp/2)\B\'^iB 



dr 
3 

2uJr, 



( 2tB'w^'^* + 2t\Bfw^'^ + 2\Bfv^''> - B'v^'^* + 2\B\-v 



2„(1) 



+ -fsm\qp/2)giB - 12T]osm\qp/2)\B\'^iB 



^p - 

dB dw^^^ 

2ujp V 
dB dw^^^ 



+ jsm\qp/2)gB - 12r]osm\qp/2)\B\''B 



2uJr 



-2tB\^'^* + 2t\B\%^'^ - A\B\V^ - B'w^'^*^ 



dB 



+ 



d^B 



dr ' 27sm2(gp/2) d^ 
3f« 



—I 



2a;p7sin {qp/2) 



+ 7 sin^ {qp/2)gB-l2r,^ sin^ {qp/2)\B\^B 

81 



4|5p— + B^ 



d^ 



d^ J ^ujl^sm\qp/2) 



\B\'B 



. (5.11b) 



In matrix form, these equations become 



+i 



dB 



+ 



d^B 



dr ' 27sin2(gp/2) d^ 



V J 

2ujp-i^ixi\qp/2) V di ) ^uh^in\qp/2) 



+ ^sm\qp/2)gB - l2riosm\qp/2)\B\''B 



81 



B\^B 



(5.12) 



The solvability condition of Eq. ()5.12|) determines the amplitude equation for B. Again, 



the Fredholm alternative theorem states that the zero eigenvalue of the operator D, 
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must be orthogonal to the right hand side of the equation. Therefore, B must satisfy 

— = ^sin\q,/2)gB+ . / ^ - 12vosin\q,/2)\B\'B 
or 27sm (gp/2) o^-^ 

+ 9 /X 4^^— + :y——\B\^B. 5.13 

2^p7sin2(gp/2) V J Su^^^ sm\qp/2y ' ^ ' 

After applying one last set of rescaling transformations, 

i-l ^5 



32 u;2^2^ sin^" (gp/2) ' 8 ^pr^oT sin'' (gp/2) ' 

16 32 
I^P ^ ^^P^o7sin^ (9p/2) and ^ ^ y o^Jr^g sin® {qp/2)g, 



we end up with an amplitude equation governed by a single parameter, 

dB d^B 2 / , ,.dB .dB*\ , „ , , , 

This amplitude equation which captures the slow dynamics of the coupled oscillators just 
above the onset of the parametric oscillations, is our central result. It is a variant of the 
CGLE [24], supplemented with the uncommon nonlinear gradient terms. Brand et. al. [28- 
30] studied similar equations that appear in the theory of binary fluid mixtures. However in 
their work all the terms in Eq. ()5.15|1 had complex coefficients, yielding localized solutions 
on which the influence of the nonlinear gradients was studied. Here we focus on extended 
solutions and show that the nonlinear gradient terms yield the wave number detunings of 
single mode solutions that, as far as we are aware of, were previously introduced by hand. 
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Chapter 6 



Single Mode Oscillations 



6.1 Single Mode Solution of The Amplitude Equation 

Once we have obtained the amphtude equation ()5.15|1 it can be used to study a variety 
of dynamical solutions, ranging from simple single-mode to more complicated nonlinear 
extended solutions, or possibly even localized solutions. Here we focus on the regime of 
small reduced amplitude g and look upon the saturation of single-mode solutions of the 
form 

B = bke-'''^, with bk = |fefc|e*^, (6.1) 
corresponding to a standing wave, determined by Eq. ()5.3|) to be 

Uq = 4|6fc| cos{upt — 7r/4) cos(gx — n/A — ip). (6.2) 

The shifted wave number q is given by 

g = g, + ^ "^^°^""'^^^/'^ ev^fc. (6.3) 

3 \Vg\ 
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The original boundary conditions u{0,t) = u{N + l,t) = impose ip = 7r/4 and require 
that the wave numbers q be quantized according to Eq. (j3.Hl . 

We substitute ()6.H) into the amplitude equation ()5.15|) and obtain 



From the linear terms in Eq. ()6.4|) we find, as expected, that for g > k"^ the zero- 
displacement solution is unstable to small perturbations of the form of ()6.HI . defining 
the parabolic neutral stability curve, shown as a dashed line in Fig. 16.21 The nonlinear 
gradients and the cubic term take the simple form 2{k — l)\bk\'^bk. For k < 1 these terms 
immediately act to saturate the growth of the amplitude assisted by the quintic term. 
Standing waves therefore bifurcate supercritically from the zero-displacement state. For 
k > 1 the cubic terms act to increase the growth of the amplitude, and saturation is 
achieved only by the quintic term. Standing waves therefore bifurcate subcritically from 
the zero-displacement state. It should be noted that similar wave-number dependent bifur- 
cations were also predicted and observed numerically in Faraday waves [26,31,32], though 
in these works the wave number detuning k was introduced a priori and was not obtained 
directly from the amplitude equations. Here the detuning is a direct consequence of the 
nonlinear gradient terms in Eq. and is explicitly related to the physical parameters 

of the system. 

The saturated amplitude \bk\, obtained by setting Eq. ()6.4|1 to zero, is given by 



In Fig. 16.11 we plot as a function of g for two values of k. The solid (dashed) lines 
are stable (unstable) states. For A; < 1 only the positive square root branch of Eq. (j6.5|) is 
obtained, and the amplitude of oscillations increases continuously from zero for g > k"^. A 



dbk 
dr 



(g-k^) bk + 2{k - l)\bk\% - \bk\%. 



(6.4) 




{k~l)± ^/{k~ly + {g-k^). 



(6.5) 
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subcritical bifurcation is obtained for A; > 1, the negative square root branch is obtained 
in addition to the positive branch for reduced driving amphtudes ranging from fc^ down to 
the saddle-node point given by 

gsn = 2k-l. (6.6) 

Since two stable solutions are obtained in this range of reduced driving amplitudes for 
subcritical bifurcations, the system exhibits hysteresis for quasistatic driving amplitude 
sweeps. Quasistatic sweeps mean that the change in the driving amplitude is much slower 
than the time required for transients to fade so that a steady-state is obtained between 
each change of the driving amplitude. 

6.2 Comparison With The Exact Form of Single Mode 
Solutions 

We now wish to compare the exact form of single mode solutions (|3.14|1 , obtained by Lifshitz 
and Cross, with the solutions obtained in the previous section. In the limit of driving 
amplitudes just above threshold and r] <^ 1, the exact amplitudes of oscillations ()3.14b|l 
correspond to Eq. (I6.5jl . The excited mode m is the mode which minimizes hm,{ujp), and its 
wave number and frequency Um are identified (up to 0(-^) corrections) with qp and ujp 
respectively. The frequency detuning fi^ corresponds to e\/5vgk, with ilc equivalent to k = 
1. This implies that for a truly continuous extended system the standing waves will always 
bifurcate supercritically with a wave number qp as long as g is increased quasistatically from 
zero. It is the discreteness of the normal modes which provides the detuning parameter 
essential for a subcritical bifurcation for a quasistatic increase of the driving amplitude. In 
Fig. 16.11 (B) we compare two exact single mode solutions calculated by Eq. ()3.14bp with 
the solution of Eq. (16. 5j) for k = 1.55. It is demonstrated that the agreement between the 



37 




Figure 6.1: The amplitude of oscillations of the single mode state |6fcp plotted as a function 
of the reduced driving amplitude g for two values of k. Solid and dashed black lines are the 
positive and negative square root branches of the calculated response in (|6.5p . the latter 
clearly unstable. (A) k = 0.9 yields a supercritical bifurcation from the zero displacement 
state into the single mode state with a continuously increasing amplitude. (B) For k = 1.55 
a subcritical bifurcation occurs at g = k"^. Colored lines mark the single mode solutions 
calculated by Eq. (j3.14bj) for e = 0.01 (blue lines) and e = 0.001 (green line). Both solutions 
are calculated for 6 = 0.01 and the rest of the parameters are chosen to yield k = 1.55. A 
good agreement with Eq. (|6.5|) is verified for e <^ ^/6. 
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amplitude equation approach solution and the exact solution is better for smaller t/\f5, as 
discussed in section W?A 

For frequency detunings and amplitudes of oscillations of order ev^, the phase obtained 
from Eq. ()3.14cp is ip = it/A + O(ev^), in agreement with the phase obtained from the 
amplitude equation. 

6.3 Secondary Instabilities 

We study secondary instabilities of the single mode solutions by performing a linear stability 
analysis of the solution ()6.H) . We substitute 

B = hkt-'^^ + (/3+e-*('=+'3)« + /3* e-*^'^-^)^) , (6.7) 

with 1/5-1- 1 -C 1, into the amplitude equation ()5.15p and linearize in I3±. Since the amplitude 
equation ()5.15|l is invariant under phase transformations B fie"""^, the stability of the 
single mode solution cannot depend on the phase of hk- We therefore assume that hk is 
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real, and linearize the following terms of Eq. ()5.15j) 



w 

\B\''B 

\B\''B 
dB 



\B\ 



B 



>dB 
dB* 
,dB* 



hi + h, + /?_)e-^«« + (/?; + P*_)e^Q^) ■ 

bl + bl (2(/5+ + /5_)e-*«« + 2(/5; + /?:)e*««) ; 
e~*'« {hi + ((2/5+ + /?_)e-''«« + (/?; + 20e^'3^)) ; 
e-'^^ {bl + bl ((3/5+ + 2/5_)e"''5€ + (2/5; + 3[3*_)e'^^)) ; 
e"*'^^ (-zA;6fc - t{k + Q)P+e~"^^ - i{k - Q)P*__e'^^) ; 

e-'^^ {-ikbl - ibl {{{2k + Q)(3+ + fc/5_)e-*«« + {k^l + {2k - Q)/3:)e^««)) ; 
e'''^ {ikbk + «(A; + Q)/3;e''3€ + ^(A; - Q)/?_e-''3€) ; 

e-'""^ {ikbl + ((2A;/3+ + (A; - Q)/5_)e-'«« + ((A; + Q)/?; + 2fc/3:)e^««)) . 

(6.8) 



The terms of order one of the equation obtained from the linearization of Eq. recover 
the same Eq. (|6.4I1 for bk- The terms with spatial dependence of e~**^^ must satisfy 

dp. 



dr 



g(3+ -{k + Qfl3+ + hi (4 {{2k + Q)I3+ + k(3_) - 2kf3+ - {k - Q)f3. 
-2{2p+ + p.)bl-{3P+ + 2p.)bt 
g-{k + Qf + bl{A{k - 1) - ?>bl + 8g/3)l /?+ + 2bl\k~l-bl + g/sl /?_. 



(6.9) 
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The terms with spatial dependence of e**^^ similarly obey 



d/3. 
'd7 



-2{P+ + 2p.)bl-{2P+ + 3p.)bl 



g-{k-Qf + hl{A{k - 1) - 36^, - 8g/3) 



Using Eq. (16 .^j] we find that (for dbk/dr = 0) 



(3. + 2bl 



k-l-bl-Q/3 



(6.10) 



g-e + A{k- l)bl - 3bl = 2bl{k -1-b 



(6.11) 



and write Eqs. ()6.9|) and ()6.in|l in matrix form as 



d_ 



M 



where 



M = 



( 



\ 



2bl{k - 1 - _ g2 _ 2Q^k - 461/3) 26i(A; -l-bl + Q/3) 

2bUk -1-bl- Q/3) 2bl{k -l-bl)-Q^ + 2Q{k - Abl/3) 



(6.12) 
\ 



We express the eigenvalues of the matrix M using its trace trM and its determinant |M| 



trM 



trMV 



\M\. 



(6.13) 



The single mode solution ()6.1|) is a stable solution only if all wave numbers Q yield negative 
eigenvalues for M, obtained for trM < and \M\ > 0. A negative trace of the matrix M 
is obtained for all wave numbers Q if 



bl> k-1 and bl > 0. 



(6.14) 
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Thus the negative square root branch in ()6.5|1 obtained for subcritical bifurcations A; > 1, 
is confirmed to be always unstable. The stability of the positive square root is determined 
by the constraint on the determinant of M, 



|M| = {2bl{k-l-bl)-Q'Y-AQ'{k-Abl/3r-Abt{k-l-bir + AbtQy9 
= Q^- 4:Q%l{k -1-bl)- AQ^e - 8blk/3 + 166^/9) + 46^QV9 
= lQ'{lQ'-bt + ^-^bl-lk^^>0. (6.15) 

In order to obtain stable single modes, this inequality should be satisfied for all wave 
numbers Q > 0^. Thus the inequality ()6.15p can be replaced with 

bt - ^-^bl + \e < ^gL., (6.16) 

where 

3 1 I li 
8ujpr]o'ysm%qp/2)eV6N + 1 

is the mode to mode separation in our rescaled units. Stable single mode oscillations 
therefore must have amplitudes obeying 



bl>0 for k< 1, (6.18a) 
bl> k-1 for k> 1, (6.18b) 



bl>— -2^' + 8^Ln, and (6.18c) 



^Perturbations with the same wave number as of the single mode, corresponding to Q = 0, decay as long 
as inequality Il(i.l4|l holds. This confirms the stability of the single mode solution towards perturbations 
with the same wave number mentioned in chanter mi Since Q = is a degenerate case of the ansatz Hfi.7|l . 
we can take /3- — and consider only the sign of A; — 1 — as apparent from Eq. 116. 9|l . The single mode 
solution is therefore stable for the positive square root branch of Eq. Hfi.5|l (5^ > fc — 1), and unstable for 
the negative branch {b^ < k — 1). 
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,2 5A; + 3 f5k + 3Y 3,, 3^^ 

The right hand side of inequality (j6.18d|) bounds the amphtude of stable single mode 
oscillations from above. The lower bound of the stable amplitudes is determined by either 
one of the inequalities ()6.18a,jl - (l6.18cl) . depending on values of k and Qmin- For k < 1 the 
right hand side of inequality (|6.18c|l determines the lower bound if 



5k + 3 Ifbk + SY 3_ 3 



+ ^QLn > 

■^2 ^ g^min 

^ 1^1 > (6-19) 



while for A; > 1 the condition states 



5A; + 3 //5A; + 3\^ 3,, 3 



A J 2 

^ (fc - If - ^{5k + 3){k-l) + - ^QLn > 

3, 5\ ^'2 ^ 



^ -(A; -l)\^-k + -j+ -k' - -Qi,^ > 

5 3 /g, " ' 



■^1<A;<--- 



2 2 V 2 



'mm 



(6.20) 



In the experimental scheme we focus on (see chapter Ej) the bifurcation from the zero dis- 
placement state to single mode oscillations always occur for \k\ < Qminf^- The stability of 
the single mode oscillations is therefore bounded from below by ()6.18a|) and (|6.18b|) . Using 
Eq. we can write the boundaries of stable single mode oscillations as a function of the 
reduced driving amplitude g and k. The lower bounds determined by inequalities ()6.18a,|l 
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and (j6.18bf) then take the simple form 



g = k'^ and g = 2k - I (6.21) 

respectively. We thus obtain that for a supercritical (subcritical) bifurcation, the positive 
square root branch of Eq. is stable for reduced driving amplitudes g just above the 
neutral stability (saddle-node) point. In Fig. 16. 21 we show the stability boundaries of single 
mode oscillations in the ((?, k) plane, describing the so-called stability balloon [22] of the 
single mode state. The solid black curve is the stability balloon obtained for a particular 
set of parameters, such that Qmm/S > 1, and thus it is bounded from bellow by Eqs. (|6.2H] . 
Outside the curve the oscillations undergo instabilities with Q = Qmin- In the limit of very 
large arrays Qmin 0, thus the single-mode state is unstable towards perturbations with 
Q ^ 0, typical of Eckhaus [33] instabilities in extended systems. In this limit the stability 
balloon shrinks to the blue curve in Fig l6.2[ 



44 



-8 -6 -4 -2 2 4 6 8 

k 

Figure 6.2: Stability boundaries of the single-mode solution of the amplitude equa- 
tion in the g vs. k plane. Dashed line: neutral stability boundary g = k"^ . Solid 
lines: stability boundary of the single-mode solution (16. Ij) as obtained by Eqs. (I6.18|) for 
a continuous spectrum (blue curve) and for = 92 with the same parameters given in 
section Em (black curve). Green line: saddle node point for the subcritical bifurcation 
gsn = 2k — 1 which coincides with the discrete stability boundary for A; > 1. 
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Chapter 7 

Numerical Simulations of Possible 
Experiments 

The purpose of our numerical work is twofold: to test our analytical predictions and to 
simulate possible experiments that may verify our results. The equations of motion ()2.4|1 are 
integrated numerically using the fourth-order Runge-Kutta method. Most of our numerical 
simulations mimic an experiment in which the ac component H is increased quasistatically 
from zero, until moderate oscillations are obtained, and then swept back to zero. This is 
done by fixing the parameters F, A, ujp, r] and the number of oscillators in the array A^, 
and scanning the values of H. Starting with driving amplitudes lower than the threshold 
2e7Co'p/A^, the initial conditions are taken to be those of the zero-displacement state Un(t = 
0) = Unit = 0) = for all n, superimposed with small random noise. The displacements of 
the beams Unit) develop in time according to the equations of motion (|2.4|1 . Once a steady 
state is obtained, the time average of the squared amplitude of each beam is calculated. 
The driving amplitude H is then increased, and the equations are integrated for the new 
H, using the displacements and velocities of the beams just before H was increased as the 
new initial conditions. Again a small random noise is added to the initial conditions in 
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order to prevent the system from staying on an unstable solution of the equations, solutions 
which are not accessible in experiment. By repeating this procedure for increasing drives, 
the amplitude of oscillations as a function of the driving amplitude is obtained, and by 
switching to decreasing sweeps of H after the array is excited, the hysteretic character of 
subcritical bifurcations can be observed. 

We test numerically our two main predictions of the amplitude equation ()5.15p : the 
existence of a wave-number dependent bifurcation and the stability boundaries of single 
mode oscillations. 



7.1 Wave-Number Dependent Bifurcation 

The first mode to emerge when increasing the driving amplitude from if = is the mode 
which minimizes Eq. ()3.12j) . whose wave number is the closest to among the spectrum 
of vibrational modes. In order to obtain the infiuence of the wave number detuning k on 
the type of bifurcation from a zero-displacement state to single-mode oscillations, the 
experimenter must have a way to control q^- By changing the number of oscillators 
in the array, the spectrum of vibrational modes can be modified, yielding different wave 
number detunings for the same driving frequency Up. Since and can differ by up to 
7i/{N + 1), \k\ is bounded from above by 

\k\ < (7.1) 

where Qmm is given by Eq. ()6.17|1 . In the limit of vary large arrays N ^ oo, q^ qp and 

the bifurcations are always supercritical^ Thus if one wants to prevent hysteresis effects in 

an experiment, large arrays should be considered. The weaker the dissipation of the system, 

the larger the array should be. For smaller arrays such that Qmin ^ 2, either supercritical 
^More accurately Qm ^ Qp + O(e^), as explained in section HaI 
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bifurcations {k < 1) or subcritical bifurcations {k > 1) can be obtained, depending on the 
specific values of Up and A. 

Our analytical treatment is valid for e, 5 -C 1 and ^ 1. We therefore choose to 
integrate arrays of about 100 oscillators with dissipation coefficients T = 0.01 and = 0.1, 
corresponding to e = 0.01 and 6 = 0.01 respectively. Setting the dc component to = 
0.25, ensures that the spectrum of vibrational modes satisfy Qmin — 2, and thus both types 
of bifurcations can be obtained, depending on the number of oscillators in the array. We 
choose the driving frequency 



to obtain = Qp for an array of = 100 oscillators and a mode number m = 73. 

We numerically integrate the equations of motion ()2.4|1 with the above parameters for 
arrays with = 100, = 98 and A^ = 92, corresponding to wave number detunings of 
/c = 0, A; ~ —0.81 and k ^ 1.55 respectively. For each array the scenario of increasing 
sweeps of the drive H followed by decreasing sweeps as described above is performed. For 
each driving amplitude, the Fourier components of the steady-state solution are computed 
to verify that only single modes are found, suggesting that in this regime of parameters 
only these states are stable. In Fig. 17.11 we plot |6fcp as a function of the reduced driving 
amplitude g for the three wave number shifts k. The symbols are the numerically computed 
Fourier components, blue x marks are calculated for increasing driving amplitude sweeps 
and red circles are calculated for decreasing sweeps, showing clear hysteresis for A; > 1. 
The solid (dashed) lines are the stable (unstable) solutions of Eq. ()6.5p . The disagreement 
between the analytic curves and the numeric integration lies within the scope of accuracy 
that the amplitude equations approach as discussed previously {g for instance is determined 




(7.2) 



up toO(^)~0.1). 
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In experiment it might be easier to control k by changing the dc component of the 
potential difference between the beams, thus changing qp rather than g^- In Fig. l7.2l we plot 
the response of an array of = 100 oscillators for two different dc components, A = 0.491 
and A = 0.5015. The values of all other parameters are as given above. Changing the 
dc component of the potential difference might however change other parameters of the 
array, such as the dissipation coefficients (see chapter HI) or the coefficient of the cubic 
elastic restoring force [8]. In addition, a change of qp involves a change in the scalings 
performed in our analytic treatment. Thus, for the sake of simplicity we have concentrated 
on controlling k by changing the number of oscillators N, but in experiment this will be 
difficult to do. 



7.2 Secondary Bifurcations 

The linear stability analysis of the standing wave state performed in section 16.31 predicts 
a transition to a new standing wave with a wave-number shift of a multiple of Tr/{N + 1) 
once the driving amplitude is increased and has crossed the upper bound of the stability 
balloon. Since the upper bound monotonically increases with k, the new wave number 
will always be larger. A sequence of three transitions, obtained numerically, is shown in 
Fig. 17. 3| superimposed with our theoretical predictions. The sequence of transitions is 
also sketched for comparison within the stability balloon in Fig. 17.41 All wave-number 
shifts obtained numerically are of 7r/(A^ + 1). Secondary instabilities are therefore another 
scenario for the experimental observation of hysteresis as a function of the applied driving 
amplitude. Once a transition has occurred, the system will return to its original state only 
when reducing the driving amplitude below the saddle node point. 
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Figure 7.1: Response of the oscillator array plotted as a function of reduced amplitude g 
for three different wave number shifts determined by fixing the number of oscillators A^. 
(A) k ~ -0.81 obtained for N = 92, (B) k = obtained for N = 100, and (C) k ~ 1.55 
obtained for N = 98. Solid and dashed lines are the positive and negative square root 
branches of the calculated response in ()6.5|1 . the latter clearly unstable. The symbols are 
numerical values obtained by numerical integration of the equations of motion ()2.4|) . Blue 
x-marks indicate solutions obtained for increasing sweeps of g, while red circles are for 
decreasing sweeps. The curves in (A) and (B) show a supercritical bifurcation, while that 
of (C) exhibits a subcritical bifurcation with clear hysteresis. The values of the parameters 
of the array are A = 0.5, ujp = 0.767445, e7 = 0.01 and rj = 0.1. 



50 



C\J 





^ ' pQQQ®>OA)< 



2.5 



g 



Figure 7.2: Response of an array of = 100 oscillators plotted as a function of reduced 
amplitude g for two wave-number shifts determined by changing the dc component of the 
potential difference between the beams. (A) k ~ —0.84 obtained for A = 0.4991 bifurcates 
supercritically, and (B) k ~ 1.44 obtained for A = 0.5015 bifurcates subcritically. All 
other parameters are set to the same values given in Fig. 17.11 
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Figure 7.3: A sequence of secondary instabilities following the initial onset of single- 
mode oscillations in an array of 92 beams and the same parameters given in Fig. 17.11 
plotted as a function of the reduced driving amplitude g. Solid (dashed) lines are stable 
(unstable) solutions defined by ()6.5p . for k = —0.81, k = 2.90 and k = 6.60 corresponding 
to the first wave number to emerge and two shifts of the wave number by Qmin and 2Qmin 
respectively. Numerical integration of the equations of motion (j2.4|] for an upward sweep 
of g (blue x-marks) , followed by a downward sweep (red circles) exhibits a strong hysteresis 
and confirms the theoretical predictions for the stability of the single mode oscillations as 
illustrated in Fig. 17.41 
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Figure 7.4: The stability balloon given in Fig. 16.21 superimposed with vertical and hor- 
izontal arrows that mark the secondary instability transitions shown in Fig. 17.31 Each 
transition shifts the wave- number detuning k by Qmin- 
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Conclusions 



We derived amplitude equations ()4.inp describing the response of large arrays of nonlinear 
coupled oscillators to parametric excitation, directly from the equations of motion yielding 
exact expressions for all the coefficients. The dynamics at the onset of oscillations was 
studied by reducing these two coupled equations into a single scaled equation governed 
by a single control parameter. Single mode standing waves were found to be the initial 
states that develop just above threshold, typical of parametric excitation. The single mode 
oscillations bifurcate from the zero-displacement state either supercritically or subcritically, 
depending on the wave number of the oscillations. The wave number dependence originates 
in the nonlinear gradient terms of the amplitude equation, which were usually disregarded 
in the past by others in the analysis of parametric oscillations above threshold. We also 
examined the stability of single mode oscillations, predicting a transition of the initial 
standing wave state to a new standing wave with a larger wave number as the driving 
amplitude is increased. 

In this work we showed that interesting response of coupled nonlinear oscillators ex- 
cited parametrically can also be obtained for quasistatic driving amplitude sweeps, rather 
than frequency sweeps usually preformed in experiments. We proposed and numerically 
demonstrated two experimental schemes for observing our predictions, hoping to draw 
more attention of experimenters to the dynamics produced by driving amplitude sweeps. 

The results obtained by the numerical integration of the equations of motion agree 
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with our analysis, supporting the validity of the amplitude equation ()5.15|1 . We therefore 
believe that the amplitude equations we derived can serve as a good starting point for 
studying other possible states of the system. One particular interesting dynamical behavior 
that can be considered is that of localized modes, often observed in arrays of coupled 
nonlinear oscillators and in other nonlinear systems as well. The conditions for obtaining 
such modes and their dynamical properties could be studied by looking for localized states 
of the amplitude equations. Another interesting aspect that can be addressed using the 
amplitude equations we derived is the response of the array to fast (rather than quasistatic) 
driving amplitude sweeps, which should lead to more complicated nonlinear response as 
observed by LC in their work. 
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